plus sn10x, SI

Baf53cre ENS neurons, SI data
Steady state, CTL x3 CKO(Ahr-cko) x3

loading 60k cells,
cellranger called 24,365 cells (demo); 24,604 cells (plus)

pure neurons downstream

load dependancies

load obj

GEX.seur <- readRDS("./GEX0429.seur.preAnno.rds")
GEX.seur

## An object of class Seurat 
## 22951 features across 13503 samples within 1 assay 
## Active assay: RNA (22951 features, 1040 variable features)
##  3 dimensional reductions calculated: pca, tsne, umap

color.FB <- ggsci::pal_d3("category20c")(20)[c(1,6,11,#16,
                                               2,7,12#,17
                                               )]

color.cnt <- color.FB[c(1,4)]

DimPlot(GEX.seur, reduction = "umap", label = T, group.by = "preAnno1") +
  DimPlot(GEX.seur, reduction = "umap", label = F, group.by = "cnt", cols = color.cnt)

filtering

GEX.seur <- subset(GEX.seur, subset = preAnno1 %in% levels(GEX.seur$preAnno1)[1:16] & DoubletFinder0.1 == "Singlet")
GEX.seur

## An object of class Seurat 
## 22951 features across 11932 samples within 1 assay 
## Active assay: RNA (22951 features, 1040 variable features)
##  3 dimensional reductions calculated: pca, tsne, umap

DimPlot(GEX.seur, reduction = "umap", label = T, group.by = "preAnno1") +
  DimPlot(GEX.seur, reduction = "umap", label = F, group.by = "cnt", cols = color.cnt)

markers.new.s <- list(EMN=c("Chat","Bnc2","Tox","Ptprt",
                       "Gfra2","Oprk1","Adamtsl1", 
                       "Fbxw15","Fbxw24","Chrna7",
                       "Satb1","Cntnap5b",
                       "Gabrb1","Nxph1","Lama2","Lrrc7",
                       "Ryr3","Eda","Tac1",
                       "Kctd8","Ntrk2","Penk",
                       "Fut9","Nfatc1","Egfr","Mgll",
                       "Chrm3"
                       ),
                 IMN=c("Nos1","Kcnab1",
                       "Gfra1","Etv1",
                       "Man1a","Airn",
                       "Adcy2","Cmah","Creb5","Vip","Pde1a",
                       "Ebf1","Gpc5"
                       ),
                 IN=c("Npas3","Synpr","St18","Gal",
                      "Kcnk13",
                      "Neurod6","Moxd1","Sctr",
                      "Piezo1","Sst","Adamts9","Kcnn2",
                      "Calb2"),
                 IPAN=c("Calcb","Nmu","Adgrg6","Pcdh10",
                        "Ngfr","Galr1","Il7","Aff2",
                        "Gpr149","Cdh6","Cdh8",
                        "Clstn2","Ano2","Ntrk3",
                        "Cpne4","Vwc2l","Cdh9","Scgn",
                        "Vcan","Cck","Piezo2","Kcnh7",
                        "Rerg","Bmpr1b","Skap1","Ntng1",
                        "Tafa2","Nxph2"),
                 CONTAM=c("Actl6b","Phox2b"),
                 pDEGs=c("Grid2","Ide","Btaf1","Cdh10","Slc15a2","Ccser1","Hdac9","Col25a1","Rspo2"))

pm.All_new.s <- DotPlot(GEX.seur, features = as.vector(unlist(markers.new.s)), group.by = "preAnno2",
        cols = c("midnightblue","darkorange1")) +
  theme(axis.text.x = element_text(angle = 60, hjust = 1, vjust = 1,size = 9.15))+ scale_y_discrete(limits=rev) + labs(title="All")
pm.All_new.s

pm.CTL_new.s <- DotPlot(subset(GEX.seur,subset=cnt %in% c("CTL")), features = as.vector(unlist(markers.new.s)), group.by = "preAnno2",
        cols = c("midnightblue","darkorange1")) +
  theme(axis.text.x = element_text(angle = 60, hjust = 1, vjust = 1,size = 9.15))+ scale_y_discrete(limits=rev) + labs(title="CTL only")
pm.CTL_new.s

pm.CKO_new.s <- DotPlot(subset(GEX.seur,subset=cnt %in% c("CKO")), features = as.vector(unlist(markers.new.s)), group.by = "preAnno2",
        cols = c("midnightblue","darkorange1")) +
  theme(axis.text.x = element_text(angle = 60, hjust = 1, vjust = 1,size = 9.15))+ scale_y_discrete(limits=rev) + labs(title="CKO only")
pm.CKO_new.s

re-clustering

GEX.seur <- FindVariableFeatures(GEX.seur, selection.method = "vst", nfeatures = 1500)

# Identify the 10 most highly variable genes
#top10 <- head(VariableFeatures(GEX.seur), 10)
top20 <- head(VariableFeatures(GEX.seur), 20)

# plot variable features with and without labels
plot1 <- VariableFeaturePlot(GEX.seur)
plot2 <- LabelPoints(plot = plot1, points = top20, repel = T)
plot1 + plot2

head(VariableFeatures(GEX.seur), 200)

##   [1] "Mgat4c"        "Cntn5"         "Zfp804a"       "Klhl1"        
##   [5] "Wdr17"         "Gm30613"       "Cntnap2"       "Kcnb2"        
##   [9] "Adgrg6"        "Robo2"         "Gal"           "Brinp3"       
##  [13] "Cpne4"         "Lingo2"        "Cdh8"          "Zfp804b"      
##  [17] "Nmu"           "Nrxn3"         "Ntng1"         "Kcnip4"       
##  [21] "Gm21847"       "Gm32647"       "Gm29521"       "Ano2"         
##  [25] "Pcdh9"         "Gm15261"       "Pdzrn4"        "Grm5"         
##  [29] "Sgcz"          "Tmeff2"        "Prkg1"         "Ntrk3"        
##  [33] "Cdh9"          "Gpr149"        "AC150683.1"    "Nrg1"         
##  [37] "Ebf1"          "Cdh18"         "Schip1"        "Dgkg"         
##  [41] "Car10"         "Cmah"          "Kcnq5"         "Fgf14"        
##  [45] "March1"        "Astn2"         "4930509J09Rik" "Vwc2l"        
##  [49] "Cdh6"          "Nxph2"         "Clstn2"        "Vip"          
##  [53] "Asic2"         "Sst"           "Dcc"           "Sema5a"       
##  [57] "Gm20754"       "Unc5d"         "Grin3a"        "Spock3"       
##  [61] "Gm26612"       "Piezo2"        "4930432L08Rik" "Dpp10"        
##  [65] "Gpc5"          "Gm15680"       "1700063D05Rik" "Efr3a"        
##  [69] "Pbx3"          "Pcdh10"        "Egfem1"        "Serpini1"     
##  [73] "Csmd3"         "Gm26737"       "Fut9"          "Ccbe1"        
##  [77] "Kcnh7"         "Zfhx3"         "Gm15584"       "Itgb6"        
##  [81] "Gm43388"       "Rasgef1b"      "Gm32916"       "Dgki"         
##  [85] "Dgkb"          "4930422I22Rik" "Plxna4"        "Myo3b"        
##  [89] "Cysltr2"       "Plcl1"         "Ttc29"         "Tbx20"        
##  [93] "1700034P13Rik" "Chrna9"        "1700111E14Rik" "9530059O14Rik"
##  [97] "Tnr"           "Chsy3"         "Cpa6"          "Pde1a"        
## [101] "Pcdh11x"       "Robo1"         "Gm20713"       "Gm49938"      
## [105] "D030068K23Rik" "Penk"          "Luzp2"         "B230323A14Rik"
## [109] "Pde4d"         "Csmd1"         "Il1rapl1"      "Ssc4d"        
## [113] "Col25a1"       "Trps1"         "Casp8"         "Trhde"        
## [117] "Myl1"          "Abca9"         "Arhgap6"       "Cacna2d3"     
## [121] "Grik1"         "Ghr"           "Rab27b"        "5530401A14Rik"
## [125] "Galntl6"       "Tafa2"         "Calcb"         "Gm11342"      
## [129] "B230312C02Rik" "Skap1"         "Gm19784"       "Kcnd2"        
## [133] "4931415C17Rik" "Tac1"          "Nwd1"          "Gm12239"      
## [137] "Galr1"         "Col8a1"        "Ccdc3"         "Gm28494"      
## [141] "Clcnka"        "Hsd17b14"      "Hs6st3"        "Gna14"        
## [145] "4930471E19Rik" "Cntn4"         "Chrm3"         "Adamts9"      
## [149] "Cux2"          "Kctd8"         "Lama2"         "Gm20560"      
## [153] "Gm26713"       "Gm50255"       "Gm37267"       "Galnt18"      
## [157] "Slc4a4"        "Nos1"          "Nxph1"         "4930587E11Rik"
## [161] "4930445B16Rik" "A330008L17Rik" "Gm49127"       "Tenm2"        
## [165] "Rarb"          "Gm49678"       "Aldh1a2"       "Cdh10"        
## [169] "2310039L15Rik" "Cadm2"         "Kif26b"        "Gabrg3"       
## [173] "Gm30094"       "Zfp286os"      "Gm11339"       "Pdgfra"       
## [177] "Gm16685"       "Nox3"          "Bmpr1b"        "Apol7b"       
## [181] "Dock10"        "Adamts6"       "Gm34544"       "Fam81b"       
## [185] "Adam2"         "Gm49953"       "Cbln2"         "Mid1"         
## [189] "Ccnb1"         "Npas3"         "Cntn3"         "C79798"       
## [193] "Kctd16"        "Atp12a"        "Gm16226"       "Tafa1"        
## [197] "Synpr"         "Cep72"         "Adap2os"       "Gm38405"

GEX.seur <- ScaleData(GEX.seur, features = rownames(GEX.seur))

## Centering and scaling data matrix

PCA

# exclude MT genes  and more 

#DIG <- grep("^Tra|^Trb|^Trg|^Trd|^Tcr|^Igm|^Igh|^Igk|^Igl|Jchain|^Hsp|^Rps|^Rpl|Hbb-|Hba-|^Dnaj|^AY|^Gm|^Hist",rownames(GEX.seur),value = T)

DIG <- grep("^mt-|^Tra|^Trb|^Trg|^Trd|^Tcr|^Igm|^Igh|^Igk|^Igl|Jchain|Mzb1|Vpreb|Lars2|Jun|Fos|^Hsp|^Rps|^Rpl|Hbb-|Hba-|^Dnaj|^AY|^AC|^AA|^AI|^AW|^BC|^Gm|^Hist|Rik$|-ps",
            rownames(GEX.seur),value = T)
CC_gene <- Hmisc::capitalize(tolower(as.vector(unlist(cc.genes.updated.2019))))


VariableFeatures(GEX.seur) <- setdiff(VariableFeatures(object = GEX.seur),
                                                c(#MT_gene,
                                                  DIG, 
                                                  CC_gene) )

GEX.seur <- RunPCA(GEX.seur, features = VariableFeatures(GEX.seur))

## PC_ 1 
## Positive:  Lrrtm4, Galntl6, Tox, Ndst4, Epha5, Kcnd2, Grik1, Bnc2, Zfp521, Pcdh7 
##     Nos1, Pde3a, Kcnab1, Lrp1b, Cdc14a, Tshz2, Chrna7, Zfp536, Syt6, Specc1 
##     Synpo2, Plcxd3, Dach1, Lrrc4c, Fbxw15, Rspo2, Brinp2, Fat3, Tenm3, Fbxw24 
## Negative:  Ano2, Ntrk3, Robo2, Tmeff2, Cdh8, Cpne4, Nrxn3, Adgrg6, Gpr149, Pdzrn4 
##     Clstn2, Cntn5, Mgat4c, Spock3, Cdh6, Zfp804a, Dgkg, Grin3a, Pcdh10, Sgcz 
##     Cysltr2, Myl1, Arhgap6, Astn2, Plxna4, Cux2, Cacna2d3, Itgb6, Iqgap2, Tnr 
## PC_ 2 
## Positive:  Rbfox1, Bnc2, Ptprt, Tshz2, Tafa1, Gpc6, Grik1, Frmd4b, St6galnac3, Brinp2 
##     Fbxw15, Caln1, Fbxw24, Pcdh7, Tox, Plcxd3, Tmem132c, Cdc14a, Tcf7l2, Specc1 
##     Pld5, Unc5c, Dock2, Adamtsl1, Oprk1, Alk, Sdk2, Nfia, Tmem132cos, Ccbe1 
## Negative:  Nos1, Auts2, Egfem1, Schip1, Kcnq5, Asic2, Cadps2, Etv1, Kcnab1, Epha5 
##     Cmah, Alcam, Dgkb, Gfra1, Dach1, Kcnt2, Hs6st3, Rgs6, Ebf1, Lrrc4c 
##     Stxbp6, Cntnap5a, Il1rapl1, Rgs7, Kcnj3, Cdh11, Pde4d, Fat3, Creb5, Kitl 
## PC_ 3 
## Positive:  Cdh18, Klhl1, Kcnip4, Pbx3, Kctd8, Csmd3, Gpc6, Htr4, Cadm2, C79798 
##     Meis1, March1, Vwc2l, Sema5a, Zfhx3, Dlc1, Gabrg3, Serpini1, Zbbx, Galnt18 
##     Cdh9, Cntn3, Skap1, Csmd1, Car10, Plcl1, Khdrbs2, Tenm2, Stxbp5l, Sybu 
## Negative:  Adgrg6, Sgcd, Nos1, Nmu, Cysltr2, Gpr149, Rgs6, Grin3a, Slc4a4, Cntnap5a 
##     Itgb6, Kcnab1, Nfia, Dapk2, Ano2, Fgf13, Cpne4, Ngfr, Ccbe1, Cbln2 
##     Zfp804a, Efr3a, Dgkg, Zfp536, Gfra1, Otof, Syt2, Lhfpl2, Calcb, Pcdh10 
## PC_ 4 
## Positive:  March1, Klhl1, Sema5a, Chsy3, Vwc2l, Cdh9, Dcc, Kcnh7, Galnt18, Ebf1 
##     Ntng1, Prune2, Kcnb2, Zfhx3, Pbx3, Bcl2, Rasgef1b, Tenm4, Trhde, Cntn4 
##     Zbbx, Il1rapl1, Sdk1, C79798, Alk, Cpa6, Gal, Sez6l, Col18a1, Lncbate1 
## Negative:  Dock10, Lingo2, Kcnip4, Ndst4, Tac1, Gda, Prkg1, Sorcs1, Nxph1, Epha5 
##     Dmd, Thsd4, Kcnt2, Csmd3, Lrrtm4, Cntn5, Lrrc7, Kctd8, Lrrc4c, Lrp1b 
##     Chl1, Cntn3, Fgf13, Lin7a, Lama2, Pgm5, Hgf, Galntl6, Pld5, Ccdc60 
## PC_ 5 
## Positive:  Trhde, Cntn4, Ptprd, Trps1, Nrg1, Ebf1, Cpa6, Csmd1, Kcnd2, Gal 
##     Col18a1, Zmat4, Npas3, Egfem1, Rmst, Luzp2, Kcnk13, Fstl4, Csmd3, Nkain3 
##     Nav2, Chsy3, Synpr, Hs6st3, Cntn3, Shisa6, Adgrl2, Myo1e, Lrp1b, Hcn1 
## Negative:  Dgkb, Alcam, Kcnt2, Klhl1, Thsd7b, Vwc2l, Dpp10, Il1rapl1, Cdh9, Mgat4c 
##     Rasgef1b, Pbx3, Fgf13, Nos1, Lrrc4c, Zfhx3, Galnt18, Auts2, Khdrbs2, Slc44a5 
##     Zbbx, C79798, Vcan, Olfr78, Rgs6, Stxbp6, Nmur2, Cdh20, P3h2, Hmcn1

length(setdiff(VariableFeatures(object = GEX.seur),
                                                c(#MT_gene,
                                                  DIG,CC_gene)))

## [1] 1027

setdiff(VariableFeatures(object = GEX.seur),
                                                c(#MT_gene,
                                                  DIG,CC_gene))[1:300]

##   [1] "Mgat4c"     "Cntn5"      "Zfp804a"    "Klhl1"      "Wdr17"     
##   [6] "Cntnap2"    "Kcnb2"      "Adgrg6"     "Robo2"      "Gal"       
##  [11] "Brinp3"     "Cpne4"      "Lingo2"     "Cdh8"       "Zfp804b"   
##  [16] "Nmu"        "Nrxn3"      "Ntng1"      "Kcnip4"     "Ano2"      
##  [21] "Pcdh9"      "Pdzrn4"     "Grm5"       "Sgcz"       "Tmeff2"    
##  [26] "Prkg1"      "Ntrk3"      "Cdh9"       "Gpr149"     "Nrg1"      
##  [31] "Ebf1"       "Cdh18"      "Schip1"     "Dgkg"       "Car10"     
##  [36] "Cmah"       "Kcnq5"      "Fgf14"      "March1"     "Astn2"     
##  [41] "Vwc2l"      "Cdh6"       "Nxph2"      "Clstn2"     "Vip"       
##  [46] "Asic2"      "Sst"        "Dcc"        "Sema5a"     "Unc5d"     
##  [51] "Grin3a"     "Spock3"     "Piezo2"     "Dpp10"      "Gpc5"      
##  [56] "Efr3a"      "Pbx3"       "Pcdh10"     "Egfem1"     "Serpini1"  
##  [61] "Csmd3"      "Fut9"       "Ccbe1"      "Kcnh7"      "Zfhx3"     
##  [66] "Itgb6"      "Rasgef1b"   "Dgki"       "Dgkb"       "Plxna4"    
##  [71] "Myo3b"      "Cysltr2"    "Plcl1"      "Ttc29"      "Tbx20"     
##  [76] "Chrna9"     "Tnr"        "Chsy3"      "Cpa6"       "Pde1a"     
##  [81] "Pcdh11x"    "Robo1"      "Penk"       "Luzp2"      "Pde4d"     
##  [86] "Csmd1"      "Il1rapl1"   "Ssc4d"      "Col25a1"    "Trps1"     
##  [91] "Casp8"      "Trhde"      "Myl1"       "Abca9"      "Arhgap6"   
##  [96] "Cacna2d3"   "Grik1"      "Ghr"        "Rab27b"     "Galntl6"   
## [101] "Tafa2"      "Calcb"      "Skap1"      "Kcnd2"      "Tac1"      
## [106] "Nwd1"       "Galr1"      "Col8a1"     "Ccdc3"      "Clcnka"    
## [111] "Hsd17b14"   "Hs6st3"     "Gna14"      "Cntn4"      "Chrm3"     
## [116] "Adamts9"    "Cux2"       "Kctd8"      "Lama2"      "Galnt18"   
## [121] "Slc4a4"     "Nos1"       "Nxph1"      "Tenm2"      "Rarb"      
## [126] "Aldh1a2"    "Cdh10"      "Cadm2"      "Kif26b"     "Gabrg3"    
## [131] "Zfp286os"   "Pdgfra"     "Nox3"       "Bmpr1b"     "Apol7b"    
## [136] "Dock10"     "Adamts6"    "Fam81b"     "Adam2"      "Cbln2"     
## [141] "Mid1"       "Ccnb1"      "Npas3"      "Cntn3"      "C79798"    
## [146] "Kctd16"     "Atp12a"     "Tafa1"      "Synpr"      "Cep72"     
## [151] "Adap2os"    "Loxl1"      "Galnt13"    "Met"        "Ano5"      
## [156] "Otof"       "Ngfr"       "Alcam"      "Rbfox1"     "Flrt2"     
## [161] "Lrrc4c"     "Grm7"       "Zbbx"       "Bmp5"       "Hccs"      
## [166] "Mme"        "Scnn1a"     "Hcn1"       "Acp7"       "Iqgap2"    
## [171] "Cadps2"     "Thsd4"      "Plxdc2"     "Eya1"       "Fbxl7"     
## [176] "Thsd7b"     "Grp"        "C1qtnf7"    "Antxr2"     "St8sia2"   
## [181] "Tenm4"      "Hsd17b3"    "Foxa3"      "Fmo2"       "Slc13a1"   
## [186] "Gpc6"       "Ptprt"      "Auts2"      "Ggt5"       "Bnc2"      
## [191] "Nell1"      "Calcrl"     "Pigk"       "Necab1"     "Ptprd"     
## [196] "Agtr1b"     "Dapk2"      "Olfr78"     "Rasgrf1"    "Prlr"      
## [201] "Fam47e"     "Kng1"       "Sorcs1"     "Hormad1"    "Ptk7"      
## [206] "Wee2"       "Dach1"      "Rerg"       "Cdhr1"      "Cpne8"     
## [211] "Hs3st4"     "Scn7a"      "B3galt1"    "Col18a1"    "Bcl2"      
## [216] "Tnfrsf23"   "Tenm1"      "Tex14"      "Hnf1aos1"   "Serpine2"  
## [221] "Tex21"      "Htr4"       "Slc2a13"    "Spag16"     "Serpinb5"  
## [226] "Qrfpr"      "Gad2"       "Cavin2"     "Aff2"       "Sorcs3"    
## [231] "Naaladl2"   "Prr16"      "Agrn"       "Esr1"       "St6galnac5"
## [236] "Mgat4e"     "Nmur2"      "Nell1os"    "Chst9"      "Prune2"    
## [241] "Slc22a20"   "Lhfpl2"     "Muc16"      "Tgfb1"      "Ror1"      
## [246] "Stxbp6"     "Rgs6"       "Lrp1b"      "Meis2"      "Cemip"     
## [251] "Zbtb7c"     "Ppp3ca"     "Pbld1"      "Mlph"       "Kif28"     
## [256] "Wt1os"      "Ovol2"      "Cga"        "Vmn2r18"    "Akr1d1"    
## [261] "Arhgap25"   "Adprhl1"    "Itgb2"      "Pnoc"       "Mast4"     
## [266] "Rmst"       "Prkca"      "Etv1"       "Vcan"       "Caln1"     
## [271] "Lmcd1"      "Hs3st2"     "Rtl4"       "Tmem132d"   "Plcb1"     
## [276] "Yae1d1"     "Sctr"       "Brinp2"     "Dmd"        "Gfra1"     
## [281] "Moxd1"      "St18"       "Adcy2"      "Nkain3"     "Adamtsl5"  
## [286] "Slc25a41"   "Fhit"       "Pla2g10"    "Epha8"      "Dach2"     
## [291] "Ndst4"      "Ildr2"      "Robo3"      "Thpo"       "Adgrl2"    
## [296] "Fbxw24"     "Sybu"       "Shisa9"     "Pde7b"      "Cntnap5b"

DimPlot(GEX.seur, reduction = "pca",dims = 1:2, group.by = "orig.ident") +
  DimPlot(GEX.seur, reduction = "pca",dims = 3:4, group.by = "orig.ident")

DimPlot(GEX.seur, reduction = "pca",dims = 1:2, group.by = "FB.info", cols = color.FB) +
  DimPlot(GEX.seur, reduction = "pca",dims = 3:4, group.by = "FB.info", cols = color.FB)

DimHeatmap(GEX.seur, dims = 1:12, cells = 1500, balanced = TRUE,ncol = 4)

decide PCs to use

ElbowPlot(GEX.seur,ndims = 50)

PCs <- 1:24
GEX.seur <- FindNeighbors(GEX.seur, dims = PCs, k.param = 20)

## Computing nearest neighbor graph

## Computing SNN

GEX.seur <- FindClusters(GEX.seur, method = 'igraph' ,resolution = 2)

## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
## 
## Number of nodes: 11932
## Number of edges: 520184
## 
## Running Louvain algorithm...
## Maximum modularity in 10 random starts: 0.8052
## Number of communities: 26
## Elapsed time: 1 seconds

Run UMAP/tSNE

GEX.seur <- RunTSNE(GEX.seur, dims=PCs, complexity = 100)
GEX.seur <- RunUMAP(GEX.seur, dims=PCs, n.neighbors = 20, seed.use = 145)

## Warning: The default method for RunUMAP has changed from calling Python UMAP via reticulate to the R-native UWOT using the cosine metric
## To use Python UMAP via reticulate, set umap.method to 'umap-learn' and metric to 'correlation'
## This message will be shown once per session

## 00:18:27 UMAP embedding parameters a = 0.9922 b = 1.112

## 00:18:27 Read 11932 rows and found 24 numeric columns

## 00:18:27 Using Annoy for neighbor search, n_neighbors = 20

## 00:18:27 Building Annoy index with metric = cosine, n_trees = 50

## 0%   10   20   30   40   50   60   70   80   90   100%

## [----|----|----|----|----|----|----|----|----|----|

## **************************************************|
## 00:18:28 Writing NN index file to temp file C:\Users\Shaorui\AppData\Local\Temp\RtmpGCQCJl\file41804e884ea1
## 00:18:28 Searching Annoy index using 1 thread, search_k = 2000
## 00:18:31 Annoy recall = 100%
## 00:18:31 Commencing smooth kNN distance calibration using 1 thread with target n_neighbors = 20
## 00:18:32 Initializing from normalized Laplacian + noise (using irlba)
## 00:18:33 Commencing optimization for 200 epochs, with 355450 positive edges
## 00:18:46 Optimization finished

DimPlot(GEX.seur, reduction = "tsne", label = T) + DimPlot(GEX.seur, reduction = "umap", label = T)

DimPlot(GEX.seur, reduction = "umap", label = T, group.by = "preAnno2") +
  DimPlot(GEX.seur, reduction = "umap", label = F, group.by = "cnt", cols = color.cnt)

DimPlot(GEX.seur, reduction = "umap", label = T, group.by = "seurat_clusters") +
  DimPlot(GEX.seur, reduction = "umap", label = F, group.by = "cnt", cols = color.cnt)

FeaturePlot(GEX.seur, reduction = "umap", features = c("nFeature_RNA","nCount_RNA","percent.mt","percent.rb"))

length(as.vector(unlist(markers.new.s)))

## [1] 92

FeaturePlot(subset(GEX.seur, subset = cnt %in% c("CTL")),
              features = as.vector(unlist(markers.new.s)), cols = c("lightgrey","red"))

FeaturePlot(subset(GEX.seur, subset = cnt %in% c("CKO")),
              features = as.vector(unlist(markers.new.s)), cols = c("lightgrey","red"))

Anno

GEX.seur$sort_clusters <- factor(as.character(GEX.seur$seurat_clusters),
                                 levels = c(0,1,2,12,11,17,9,10,18,
                                            4,5,14,8,7,13,16,23,
                                            6,19,24,
                                            3,20,
                                            15,21,
                                            25,22))

table(GEX.seur@meta.data[,c("preAnno2","sort_clusters")])

##          sort_clusters
## preAnno2     0    1    2   12   11   17    9   10   18    4    5   14    8    7
##   EMN1    1197  856  706  420  295    0    5    0    0    0    0    0    0    0
##   EMN2       0    0    0    0  145  188    7    0    0    0    0    0    0    0
##   EMN3       0    0    0    0    2  137  503   41    4    0    0    0    0    0
##   EMN4       0    0    0    0    0    1   13  482    4    0    0    0    0    0
##   EMN5       0    0    0    0    1    0    1    2  308    0    0    0    0    0
##   IMN1       0    0    0    0    0    0    0    0    0  583  564  311  542  455
##   IMN2       0    0    0    0    0    0    0    0    0    0    0   97    1   90
##   IMN3       0    0    0    0    0    0    0    0    0    0    1    0    0    0
##   IMN4       0    0    0    0    0    0    0    0    0    0    0    0    0    0
##   IN1        0    0    0    0    0    0    0    0    0    0    0    0    0    0
##   IN2        0    0    0    0    0    0    0    0    0    0    0    0    0    0
##   IN3        0    0    0    0    0    0    0    0    0    0    0    0    0    0
##   IPAN1.1    0    0    0    0    0    0    0    0    0    0    0    0    0    0
##   IPAN1.2    0    0    0    0    0    0    0    0    0    0    0    0    0    0
##   IPAN2.1    0    0    0    0    0    0    0    0    0    0    0    0    0    0
##   IPAN2.2    0    0    0    0    0    0    1    0    0    0    0    0    0    0
##   IPAN3      0    0    0    0    0    0    0    0    0    0    0    0    0    0
##   IPAN4      0    0    0    0    0    0    0    0    0    0    0    0    0    0
##   MIX        0    0    0    0    0    0    0    0    0    0    0    0    0    0
##   Glia       0    0    0    0    0    0    0    0    0    0    0    0    0    0
##   SMC        0    0    0    0    0    0    0    0    0    0    0    0    0    0
##          sort_clusters
## preAnno2    13   16   23    6   19   24    3   20   15   21   25   22
##   EMN1       0    0    0    0    0   33    0    0    0    0    0    0
##   EMN2       0    0    0    0    0    0    0    0    0    0    0    0
##   EMN3       0    0    0    0    0    0    0    0    0    0    0    0
##   EMN4       0    0    0    0    0    0    0    0    0    0    0    0
##   EMN5       0    0    0    0    0    0    0    0    0    0    0    0
##   IMN1       9    0    0    0    0    0    0    0    0    0    0    0
##   IMN2     404   22    0    0    0    0    0    0    0    0    0    0
##   IMN3       0  306   42    0    0    0    0    0    0    0    0    0
##   IMN4       0    0  212    0    0    0    0    0    0    0    1    0
##   IN1        0    0    0  539    5    0    0    0    0    0    0    0
##   IN2        0    0    0    7  301    0    0    0    0    0    0    0
##   IN3        0    0    0    0    0  118    0    0    0    0    0    0
##   IPAN1.1    0    0    0    0    0    0  611    1    0    0    0    0
##   IPAN1.2    0    0    0    0    0    0   51  299    0    0    0    0
##   IPAN2.1    0    0    0    0    0    0    0    1  294    5    0    0
##   IPAN2.2    0    1    0    0    0    0    0    0   40  287    0    0
##   IPAN3      0    0    1    0    0    0    0    0    0    0   95    0
##   IPAN4      0    0    0    0    0    0    0    0    0    0    8  276
##   MIX        0    0    0    0    0    0    0    0    0    0    0    0
##   Glia       0    0    0    0    0    0    0    0    0    0    0    0
##   SMC        0    0    0    0    0    0    0    0    0    0    0    0

markers.new.s <- list(EMN=c("Chat","Bnc2","Tox","Ptprt",
                       "Gfra2","Oprk1","Adamtsl1", 
                       "Fbxw15","Fbxw24","Chrna7",
                       "Satb1","Cntnap5b",
                       "Gabrb1","Nxph1","Lama2","Lrrc7",
                       "Ryr3","Eda","Tac1",
                       "Kctd8","Ntrk2","Penk",
                       "Fut9","Nfatc1","Egfr","Mgll",
                       "Chrm3"
                       ),
                 IMN=c("Nos1","Kcnab1",
                       "Gfra1","Etv1",
                       "Man1a","Airn",
                       "Adcy2","Cmah","Creb5","Vip","Pde1a",
                       "Ebf1","Gpc5"
                       ),
                 IN=c("Npas3","Synpr","St18","Gal",
                      "Kcnk13",
                      "Neurod6","Moxd1","Sctr",
                      "Piezo1","Sst","Adamts9","Kcnn2",
                      "Calb2"),
                 IPAN=c("Calcb","Nmu","Adgrg6","Pcdh10",
                        "Ngfr","Galr1","Il7","Aff2",
                        "Gpr149","Cdh6","Cdh8",
                        "Clstn2","Ano2","Ntrk3",
                        "Cpne4","Vwc2l","Cdh9","Scgn",
                        "Vcan","Cck","Piezo2","Kcnh7",
                        "Rerg","Bmpr1b","Skap1","Ntng1",
                        "Tafa2","Nxph2"),
                 CONTAM=c("Actl6b","Phox2b"),
                 pDEGs=c("Grid2","Ide","Btaf1","Cdh10","Slc15a2","Ccser1","Hdac9","Col25a1","Rspo2"))

pm.All_new.s <- DotPlot(GEX.seur, features = as.vector(unlist(markers.new.s)), group.by = "sort_clusters",
        cols = c("midnightblue","darkorange1")) +
  theme(axis.text.x = element_text(angle = 60, hjust = 1, vjust = 1,size = 9.15))+ scale_y_discrete(limits=rev) + labs(title="All")
pm.All_new.s

# Anno1            
GEX.seur$Anno1 <- as.character(GEX.seur$seurat_clusters)

GEX.seur$Anno1[GEX.seur$Anno1 %in% c(0,1,2,12)] <- "EMN1"
GEX.seur$Anno1[GEX.seur$Anno1 %in% c(11,17)] <- "EMN2"
GEX.seur$Anno1[GEX.seur$Anno1 %in% c(9)] <- "EMN3"
GEX.seur$Anno1[GEX.seur$Anno1 %in% c(10)] <- "EMN4"
GEX.seur$Anno1[GEX.seur$Anno1 %in% c(18)] <- "EMN5"

GEX.seur$Anno1[GEX.seur$Anno1 %in% c(4,5,14,8)] <- "IMN1"
GEX.seur$Anno1[GEX.seur$Anno1 %in% c(7,13)] <- "IMN2"
GEX.seur$Anno1[GEX.seur$Anno1 %in% c(16)] <- "IMN3"
GEX.seur$Anno1[GEX.seur$Anno1 %in% c(23)] <- "IMN4"

GEX.seur$Anno1[GEX.seur$Anno1 %in% c(6)] <- "IN1"
GEX.seur$Anno1[GEX.seur$Anno1 %in% c(19)] <- "IN2"
GEX.seur$Anno1[GEX.seur$Anno1 %in% c(24)] <- "IN3"

GEX.seur$Anno1[GEX.seur$Anno1 %in% c(3,20)] <- "IPAN1"
GEX.seur$Anno1[GEX.seur$Anno1 %in% c(15,21)] <- "IPAN2"
GEX.seur$Anno1[GEX.seur$Anno1 %in% c(25)] <- "IPAN3"
GEX.seur$Anno1[GEX.seur$Anno1 %in% c(22)] <- "IPAN4"


GEX.seur$Anno1 <- factor(GEX.seur$Anno1,
                            levels = c(paste0("EMN",1:5),
                                       paste0("IMN",1:4),
                                       paste0("IN",1:3),
                                       paste0("IPAN",1:4)
                                       ))


# Anno2            
GEX.seur$Anno2 <- as.character(GEX.seur$seurat_clusters)

GEX.seur$Anno2[GEX.seur$Anno2 %in% c(0,1,2,12)] <- "EMN1"
GEX.seur$Anno2[GEX.seur$Anno2 %in% c(11,17)] <- "EMN2"
GEX.seur$Anno2[GEX.seur$Anno2 %in% c(9)] <- "EMN3"
GEX.seur$Anno2[GEX.seur$Anno2 %in% c(10)] <- "EMN4"
GEX.seur$Anno2[GEX.seur$Anno2 %in% c(18)] <- "EMN5"

GEX.seur$Anno2[GEX.seur$Anno2 %in% c(4,5,14,8)] <- "IMN1"
GEX.seur$Anno2[GEX.seur$Anno2 %in% c(7,13)] <- "IMN2"
GEX.seur$Anno2[GEX.seur$Anno2 %in% c(16)] <- "IMN3"
GEX.seur$Anno2[GEX.seur$Anno2 %in% c(23)] <- "IMN4"

GEX.seur$Anno2[GEX.seur$Anno2 %in% c(6)] <- "IN1"
GEX.seur$Anno2[GEX.seur$Anno2 %in% c(19)] <- "IN2"
GEX.seur$Anno2[GEX.seur$Anno2 %in% c(24)] <- "IN3"

GEX.seur$Anno2[GEX.seur$Anno2 %in% c(3)] <- "IPAN1.1"
GEX.seur$Anno2[GEX.seur$Anno2 %in% c(20)] <- "IPAN1.2"
GEX.seur$Anno2[GEX.seur$Anno2 %in% c(15)] <- "IPAN2.1"
GEX.seur$Anno2[GEX.seur$Anno2 %in% c(21)] <- "IPAN2.2"
GEX.seur$Anno2[GEX.seur$Anno2 %in% c(25)] <- "IPAN3"
GEX.seur$Anno2[GEX.seur$Anno2 %in% c(22)] <- "IPAN4"


GEX.seur$Anno2 <- factor(GEX.seur$Anno2,
                            levels = c(paste0("EMN",1:5),
                                       paste0("IMN",1:4),
                                       paste0("IN",1:3),
                                       paste0("IPAN",c(1.1,1.2,2.1,2.2,3,4))
                                       ))

DimPlot(GEX.seur, reduction = "umap", label = T, group.by = "Anno1") +
  DimPlot(GEX.seur, reduction = "umap", label = T, group.by = "Anno2")

markers.new.s <- list(EMN=c("Chat","Bnc2","Tox","Ptprt",
                       "Gfra2","Oprk1","Adamtsl1", 
                       "Fbxw15","Fbxw24","Chrna7",
                       "Satb1","Cntnap5b",
                       "Gabrb1","Nxph1","Lama2","Lrrc7",
                       "Ryr3","Eda","Tac1",
                       "Kctd8","Ntrk2","Penk",
                       "Fut9","Nfatc1","Egfr","Mgll",
                       "Chrm3"
                       ),
                 IMN=c("Nos1","Kcnab1",
                       "Gfra1","Etv1",
                       "Man1a","Airn",
                       "Adcy2","Cmah","Creb5","Vip","Pde1a",
                       "Ebf1","Gpc5"
                       ),
                 IN=c("Npas3","Synpr","St18","Gal",
                      "Kcnk13",
                      "Neurod6","Moxd1","Sctr",
                      "Piezo1","Sst","Adamts9","Kcnn2",
                      "Calb2"),
                 IPAN=c("Calcb","Nmu","Adgrg6","Pcdh10",
                        "Ngfr","Galr1","Il7","Aff2",
                        "Gpr149","Cdh6","Cdh8",
                        "Clstn2","Ano2","Ntrk3",
                        "Cpne4","Vwc2l","Cdh9","Scgn",
                        "Vcan","Cck","Piezo2","Kcnh7",
                        "Rerg","Bmpr1b","Skap1","Ntng1",
                        "Tafa2","Nxph2"),
                 CONTAM=c("Actl6b","Phox2b"),
                 pDEGs=c("Grid2","Ide","Btaf1","Cdh10","Slc15a2","Ccser1","Hdac9","Col25a1","Rspo2"))

pm.A1_new.s <- DotPlot(GEX.seur, features = as.vector(unlist(markers.new.s)), group.by = "Anno1",
        cols = c("midnightblue","darkorange1")) +
  theme(axis.text.x = element_text(angle = 60, hjust = 1, vjust = 1,size = 9.15))+ scale_y_discrete(limits=rev) + labs(title="All Anno1")
pm.A1_new.s

pm.A2_new.s <- DotPlot(GEX.seur, features = as.vector(unlist(markers.new.s)), group.by = "Anno2",
        cols = c("midnightblue","darkorange1")) +
  theme(axis.text.x = element_text(angle = 60, hjust = 1, vjust = 1,size = 9.15))+ scale_y_discrete(limits=rev) + labs(title="All Anno2")
pm.A2_new.s

I have noticed that annotations of EMN1-5 each time are a little different from others ..
the boudaries are not that distinct to separate them,
better to do integration if have to measure the subtypes for multiple datasets

cell composition

GEX.seur$FB.info <- factor(as.character(GEX.seur$FB.info),
                           levels = c(paste0("CTL.",1:3),
                                      paste0("CKO.",1:3)))

GEX.seur$rep <- paste0("rep",
                        gsub("CTL|CKO","",as.character(GEX.seur$FB.info)))

DimPlot(GEX.seur, reduction = "umap", label = T, group.by = "Anno1") +
  DimPlot(GEX.seur, reduction = "umap", label = F, group.by = "cnt", cols = color.cnt)

DimPlot(GEX.seur, reduction = "umap", label = T, group.by = "Anno2") +
  DimPlot(GEX.seur, reduction = "umap", label = F, group.by = "cnt", cols = color.cnt)

heatmap.1

cowplot::plot_grid(
pheatmap::pheatmap(table(FB.info=GEX.seur$FB.info,
      clusters=GEX.seur$Anno1),
                   main = "Cell Count",
                   gaps_row = c(3),
      gaps_col = c(5,9,12),
                   cluster_rows = F, cluster_cols = F, display_numbers = T, number_format = "%.0f", legend = F, silent = T)$gtable,

pheatmap::pheatmap(100*rowRatio(table(FB.info=GEX.seur$FB.info,
                                clusters=GEX.seur$Anno1)),
                   main = "Cell Ratio",
                   gaps_row = c(3),
      gaps_col = c(5,9,12),
                   cluster_rows = F, cluster_cols = F, display_numbers = T, number_format = "%.1f", legend = F, silent =T)$gtable,
ncol = 1)

barplot.1

stat_Anno1 <- GEX.seur@meta.data[,c("cnt","rep","Anno1")]

stat_Anno1.s <- stat_Anno1 %>%
  group_by(cnt,rep,Anno1) %>%
  summarise(count=n()) %>%
  mutate(prop= count/sum(count)) %>%
  ungroup

## `summarise()` has grouped output by 'cnt', 'rep'. You can override using the
## `.groups` argument.

pcom.Anno1 <- stat_Anno1.s %>%
  ggplot(aes(x = Anno1, y = 100*prop, color=cnt)) +
  geom_bar(stat="summary", fun="mean", position = position_dodge(0.75), width = 0.58, fill="white") +
  theme_classic(base_size = 15) +
  scale_color_manual(values = c("skyblue","pink"), name="") +
  labs(title="Cell Composition of SI data, Anno1", y = "Proportion") +
  theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1, size = 9.6)) +
  
  geom_point(aes(x=Anno1, y = 100*prop, color= cnt),
             position = position_dodge(0.75),
             shape=16,alpha=0.75,size=2.15,
             stroke=0.15, show.legend=F)

pcom.Anno1

sig.test

glm.nb - anova.LRT

Sort.N <- c("CTL","CKO")

glm.nb_anova.Anno1 <- list()

for(x1 in 1:2){
  for(x2 in 1:2){
    if(x2>x1){
      
      stat_Anno1.s_N <- stat_Anno1.s %>% filter(cnt %in% c(Sort.N[x1],Sort.N[x2]))
      
      total.N <- stat_Anno1.s_N %>%
        group_by(cnt, rep) %>%
        summarise(total=sum(count)) %>% ungroup() %>% as.data.frame()
      
      rownames(total.N) <- paste0(as.character(total.N$cnt),
                                  "_",
                                  as.character(total.N$rep))
      
      stat_Anno1.s_N$total <- total.N[paste0(stat_Anno1.s_N$cnt,
                                            "_",
                                            stat_Anno1.s_N$rep),"total"]
      
      glm.nb_anova.Anno1[[paste0(Sort.N[x1],"vs",Sort.N[x2])]] <- list()
      
      for(AA in levels(stat_Anno1.s_N$Anno1)){
        glm.nb_anova.Anno1[[paste0(Sort.N[x1],"vs",Sort.N[x2])]][[AA]] <- tryCatch({
          aaa <- stat_Anno1.s_N %>% filter(Anno1==AA)
          aaa.glmnb <- MASS::glm.nb(count ~ cnt + offset(log(total)), data = aaa, maxit=1000)
          aaa.anova <- anova(aaa.glmnb, test = "LRT")
          aaa.anova$'Pr(>Chi)'[2]
        }, error=function(e){
        
          tryCatch({
            aaa <- stat_Anno1.s_N %>% filter(Anno1==AA)
            aaa.glmnb <- MASS::glm.nb(count ~ cnt + offset(log(total)), data = aaa, maxit=100)
            aaa.anova <- anova(aaa.glmnb, test = "LRT")
            aaa.anova$'Pr(>Chi)'[2]
          }, error=function(e){
            
            tryCatch({
              aaa <- stat_Anno1.s_N %>% filter(Anno1==AA)
              aaa.glmnb <- MASS::glm.nb(count ~ cnt + offset(log(total)), data = aaa, maxit=10)
              aaa.anova <- anova(aaa.glmnb, test = "LRT")
              aaa.anova$'Pr(>Chi)'[2]
            }, error = function(e){
              NA
            })     
          })
        })
      }
      glm.nb_anova.Anno1[[paste0(Sort.N[x1],"vs",Sort.N[x2])]] <- unlist(glm.nb_anova.Anno1[[paste0(Sort.N[x1],"vs",Sort.N[x2])]]) 
    }
  }
}

glm.nb_anova.Anno1_df <- t(data.frame(Reduce(function(x,y){rbind(x,y)},glm.nb_anova.Anno1)))
rownames(glm.nb_anova.Anno1_df) <- names(glm.nb_anova.Anno1)
colnames(glm.nb_anova.Anno1_df) <- gsub("X","C",colnames(glm.nb_anova.Anno1_df))
glm.nb_anova.Anno1_df

##               EMN1      EMN2      EMN3        EMN4      EMN5      IMN1
## CTLvsCKO 0.3141664 0.3103288 0.1787928 0.001014898 0.5353266 0.1853838
##               IMN2      IMN3       IMN4          IN1       IN2         IN3
## CTLvsCKO 0.1070702 0.6577923 0.05522512 4.638542e-05 0.2279499 0.009895313
##              IPAN1     IPAN2     IPAN3     IPAN4
## CTLvsCKO 0.4831668 0.8723714 0.2798343 0.8738828

round(glm.nb_anova.Anno1_df,5)

##             EMN1    EMN2    EMN3    EMN4    EMN5    IMN1    IMN2    IMN3
## CTLvsCKO 0.31417 0.31033 0.17879 0.00101 0.53533 0.18538 0.10707 0.65779
##             IMN4   IN1     IN2    IN3   IPAN1   IPAN2   IPAN3   IPAN4
## CTLvsCKO 0.05523 5e-05 0.22795 0.0099 0.48317 0.87237 0.27983 0.87388

heatmap.2

cowplot::plot_grid(
pheatmap::pheatmap(table(FB.info=GEX.seur$FB.info,
      clusters=GEX.seur$Anno2),
                   main = "Cell Count",
                   gaps_row = c(3),
      gaps_col = c(5,9,12),
                   cluster_rows = F, cluster_cols = F, display_numbers = T, number_format = "%.0f", legend = F, silent = T)$gtable,

pheatmap::pheatmap(100*rowRatio(table(FB.info=GEX.seur$FB.info,
                                clusters=GEX.seur$Anno2)),
                   main = "Cell Ratio",
                   gaps_row = c(3),
      gaps_col = c(5,9,12),
                   cluster_rows = F, cluster_cols = F, display_numbers = T, number_format = "%.1f", legend = F, silent =T)$gtable,
ncol = 1)

barplot.2

stat_Anno2 <- GEX.seur@meta.data[,c("cnt","rep","Anno2")]

stat_Anno2.s <- stat_Anno2 %>%
  group_by(cnt,rep,Anno2) %>%
  summarise(count=n()) %>%
  mutate(prop= count/sum(count)) %>%
  ungroup

## `summarise()` has grouped output by 'cnt', 'rep'. You can override using the
## `.groups` argument.

pcom.Anno2 <- stat_Anno2.s %>%
  ggplot(aes(x = Anno2, y = 100*prop, color=cnt)) +
  geom_bar(stat="summary", fun="mean", position = position_dodge(0.75), width = 0.58, fill="white") +
  theme_classic(base_size = 15) +
  scale_color_manual(values = c("skyblue","pink"), name="") +
  labs(title="Cell Composition of SI data, Anno2", y = "Proportion") +
  theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1, size = 9.6)) +
  
  geom_point(aes(x=Anno2, y = 100*prop, color= cnt),
             position = position_dodge(0.75),
             shape=16,alpha=0.75,size=2.15,
             stroke=0.15, show.legend=F)

pcom.Anno2

sig.test

glm.nb - anova.LRT

Sort.N <- c("CTL","CKO")

glm.nb_anova.Anno2 <- list()

for(x1 in 1:2){
  for(x2 in 1:2){
    if(x2>x1){
      
      stat_Anno2.s_N <- stat_Anno2.s %>% filter(cnt %in% c(Sort.N[x1],Sort.N[x2]))
      
      total.N <- stat_Anno2.s_N %>%
        group_by(cnt, rep) %>%
        summarise(total=sum(count)) %>% ungroup() %>% as.data.frame()
      
      rownames(total.N) <- paste0(as.character(total.N$cnt),
                                  "_",
                                  as.character(total.N$rep))
      
      stat_Anno2.s_N$total <- total.N[paste0(stat_Anno2.s_N$cnt,
                                            "_",
                                            stat_Anno2.s_N$rep),"total"]
      
      glm.nb_anova.Anno2[[paste0(Sort.N[x1],"vs",Sort.N[x2])]] <- list()
      
      for(AA in levels(stat_Anno2.s_N$Anno2)){
        glm.nb_anova.Anno2[[paste0(Sort.N[x1],"vs",Sort.N[x2])]][[AA]] <- tryCatch({
          aaa <- stat_Anno2.s_N %>% filter(Anno2==AA)
          aaa.glmnb <- MASS::glm.nb(count ~ cnt + offset(log(total)), data = aaa, maxit=1000)
          aaa.anova <- anova(aaa.glmnb, test = "LRT")
          aaa.anova$'Pr(>Chi)'[2]
        }, error=function(e){
        
          tryCatch({
            aaa <- stat_Anno2.s_N %>% filter(Anno2==AA)
            aaa.glmnb <- MASS::glm.nb(count ~ cnt + offset(log(total)), data = aaa, maxit=100)
            aaa.anova <- anova(aaa.glmnb, test = "LRT")
            aaa.anova$'Pr(>Chi)'[2]
          }, error=function(e){
            
            tryCatch({
              aaa <- stat_Anno2.s_N %>% filter(Anno2==AA)
              aaa.glmnb <- MASS::glm.nb(count ~ cnt + offset(log(total)), data = aaa, maxit=10)
              aaa.anova <- anova(aaa.glmnb, test = "LRT")
              aaa.anova$'Pr(>Chi)'[2]
            }, error = function(e){
              NA
            })     
          })
        })
      }
      glm.nb_anova.Anno2[[paste0(Sort.N[x1],"vs",Sort.N[x2])]] <- unlist(glm.nb_anova.Anno2[[paste0(Sort.N[x1],"vs",Sort.N[x2])]]) 
    }
  }
}

glm.nb_anova.Anno2_df <- t(data.frame(Reduce(function(x,y){rbind(x,y)},glm.nb_anova.Anno2)))
rownames(glm.nb_anova.Anno2_df) <- names(glm.nb_anova.Anno2)
colnames(glm.nb_anova.Anno2_df) <- gsub("X","C",colnames(glm.nb_anova.Anno2_df))
glm.nb_anova.Anno2_df

##               EMN1      EMN2      EMN3        EMN4      EMN5      IMN1
## CTLvsCKO 0.3141664 0.3103288 0.1787928 0.001014898 0.5353266 0.1853838
##               IMN2      IMN3       IMN4          IN1       IN2         IN3
## CTLvsCKO 0.1070702 0.6577923 0.05522512 4.638542e-05 0.2279499 0.009895313
##            IPAN1.1   IPAN1.2   IPAN2.1   IPAN2.2     IPAN3     IPAN4
## CTLvsCKO 0.5696465 0.3524691 0.8613843 0.6981254 0.2798343 0.8738828

round(glm.nb_anova.Anno2_df,5)

##             EMN1    EMN2    EMN3    EMN4    EMN5    IMN1    IMN2    IMN3
## CTLvsCKO 0.31417 0.31033 0.17879 0.00101 0.53533 0.18538 0.10707 0.65779
##             IMN4   IN1     IN2    IN3 IPAN1.1 IPAN1.2 IPAN2.1 IPAN2.2   IPAN3
## CTLvsCKO 0.05523 5e-05 0.22795 0.0099 0.56965 0.35247 0.86138 0.69813 0.27983
##            IPAN4
## CTLvsCKO 0.87388

#saveRDS(GEX.seur,"./GEX0430.seur.pure_Anno.rds")